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ABSTRACT 


This  study  examines  mixing  characteristics  of  double-diffusive  convection  for  a  wide 
range  of  fluids.  Our  approach  involves  Direct  Numerical  Simulation  (DNS)  utilizing  de- 
aliased  pseudo-spectral  method.  To  expedite  these  simulations  the  numerical  algorithm 
was  parallelized  using  Message  Passing  Interface  (MPI)  calculations  in  both  two  and 
three  dimensions.  A  theoretical  model  of  equilibrium  double-diffusive  transport  is 
presented,  which  emphasizes  the  role  of  secondary  instabilities  of  salt  fingers  in 
saturation  of  their  linear  growth.  Theory  assumes  that  the  fully  developed  equilibrium 
state  is  characterized  by  the  comparable  growth  rates  of  primary  and  secondary 
instabilities.  This  assumption  makes  it  possible  to  fonnulate  a  simple  and  efficient 
algorithm  for  computing  diffusivities  of  heat  and  salt  as  a  function  of  the  background 
property  gradients  and  molecular  parameters.  The  model  predicts  that  the  double- 
diffusive  transport  of  heat  and  salt  rapidly  intensifies  with  decreasing  density  ratio. 
Fluxes  are  less  sensitive  to  molecular  characteristics,  mildly  increasing  with  Prandtl 
number  (Pr)  and  decreasing  with  diffusivity  ratio  (x).  Theory  is  successfully  tested  by  a 
series  of  direct  numerical  simulations  which  span  a  wide  range  of  Pr  and  x. 

Double  diffusion  occurs  on  the  micro-scale  and  computer  technology  is  just  now 
reaching  the  processing  speeds  needed  to  fully  resolve  this  complex  phenomenon  in  three 
dimensions.  In  addition  to  the  well-known  “salt  fingering”  within  the  oceans,  double 
diffusion  occurs  in  many  other  statically  stable  regions,  both  terrestrially  and  beyond  our 
atmosphere.  Understanding  this  phenomenon  could  prove  essential  as  we  continue  to 
discover  new  worlds  and  new  areas  within  our  galaxy,  including  here  on  terra  finna. 
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I.  INTRODUCTION 


Double  diffusion  is  the  instability  of  a  stratified  fluid  at  rest  whose  density  is 
determined  by  two  components  diffusing  at  different  rates.  Such  a  configuration  can  be 
unstable  even  if  the  density  of  the  fluid  is  increasing  downwards.  The  resulting  double- 
diffusive  convection  has  long  been  recognized  as  a  significant,  and  in  many  cases 
dominant,  mixing  process  in  the  ocean,  where  wann  and  salty  water  is  often  located 
above  cold  and  fresh.  In  this  case,  the  faster  diffuser  (temperature  T)  is  stabilizing  and  the 
slower  diffuser  (salinity  S)  is  destabilizing,  resulting  in  the  salt  finger  form  of  double- 
diffusive  convection  (Figure  1),  which  is  the  main  focus  of  this  discussion. 
Equation  Section  (Next) 

Salt  fingers  are  mere  millimeters  in  diameter  and  can  grow  to  a  length  of  over  20 
centimeters  (Pickard  1990).  Double-diffusive  convection  interactions  were  originally 
noted  by  W.  S.  Jevons  (Jevons  1857)  and  were  the  first  recorded  occurrences  of  this 
phenomenon.  At  the  time,  Jevons  was  trying  to  model  cirrus  clouds  when  he 
inadvertently  discovered  double-diffusive  convection  within  the  laboratory  environment. 
In  his  entry  to  Edinburgh  and  Dublin  Philosophical  Magazine  and  Journal  of  Science  in 
1857,  Jevons  used  a  warm  sugar  solution  over  a  fresh  water  solution  in  a  lab  experiment 
and  reported  below: 

The  parts  of  these  strata,  however,  which  are  immediately  in  contact,  soon 
communicate  their  heat  and  tend  to  assume  a  mean  temperature;  and  it  is 
evident  that  whenever  this  is  the  case,  the  portions  of  liquid  containing 
sugar  must  always  be  slightly  denser  than  those  that  are  pure,  and  must 
consequently  sink  below  and  displace  the  latter. 

We  shall  thus  have  portions  of  the  upper  stratum  continually  sinking  into 
the  lower,  and  corresponding  portions  of  the  lower  rising  through  the 
upper;  and  this  movement,  as  the  experiment  demonstrates,  takes  place  by 
an  interfiltration  of  minute,  thread-like  streams. 

The  first  physical  model  of  salt  fingers  was  developed  by  Dr.  Melvin  E.  Stern 
(Stern  1960)  more  than  a  century  later,  which  marks  the  beginning  of  the  modern  theory 
of  double-diffusive  convection.  This  phenomenon  is  better  understood  today  due  to  the 

use  of  high  performance  computers  (Radko  and  Stern  1999;  Kimura  and  Smyth  2007; 
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Prikasky  2007;  Caro  2009)  important  insights  were  brought  by  laboratory  studies 
(Krishnamurti  2003)  and  field  measurements  with  in  the  ocean  (Gargett  1982).  While  it  is 
still  not  fully  understood  how  these  minute  phenomena  impact  larger  scale  circulations,  it 
is  widely  believed  that  they  may  hold  the  key  to  understanding  the  driving  forces  behind 
the  general  circulation  of  the  ocean  and  therefore,  the  Earth’s  climate. 

In  addition  to  oceanographic  applications,  compositionally  driven  double- 
diffusive  convection  affects  heat  and  material  transport  in  a  variety  of  other  geo-  and 
astrophysical  fluid  systems,  from  magmatic  melts  (Tait  and  Jaupart  1989)  to  the  interiors 
of  giant  planets  and  stars  (Guillot  1999;  Vauclair  2004;  Charbonnel  and  Zahn  2007; 
Stancliffe  et  al.„  2007). 

A.  SALT-FINGERS  IN  THE  OCEAN 

Salt  fingers  occur  in  nature  in  statically  stable  environments  (Figures  2  and  3), 
where  temperature  and  salinity  perturbations  diffuse  at  different  rates  causing  fingers  to 
form  within  a  body  of  water.  Salt  fingers  are  common  in  the  main  thermocline,  where 
water  tends  to  be  wanner  and  saltier  near  the  sea-surface  due  to  a  combination  of  surface 
evaporation  and  heating.  Individual  parcels  displaced  downward  lose  their  thermal 
properties  faster  than  salinity,  causing  them  to  be  denser  than  the  surrounding  fluid 
(Figure  1)  and  continue  to  sink.  This  causes  the  parcel  to  continue  descending  (Stern 
1960;  Schmitt  et  al.„  2005),  which  results  in  vigorous  double-diffusive  convection 
discussed  in  this  study. 
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Salt  Fingers 


Figure  1.  Schematic  diagram  illustrating  the  physical  mechanisms  of  double-diffusive 

convection. 

An  example  of  this  process  can  be  found  at  the  output  of  the  Mediterranean 
through  the  Straits  of  Gibraltar.  The  outflow  from  the  Mediterranean  Sea  is  much  saltier 
and  wanner  than  the  Atlantic  Ocean  water  at  depth  (approximately  2400  meters)  causing 
double-diffusive  interactions  deep  underneath  the  surface  of  the  ocean. 

It  is  also  believed  that  salt  fingers  contribute  to  the  ocean’s  large  scale  mixing  and 
could  be  important  to  climate  models  (Gargett  and  Schmitt  1982).  Today’s  climate 
models  cannot  resolve  such  minute  scales  within  a  reasonable  forecast  time,  so  finding 
physical  parameterizations  that  can  be  fed  into  these  models  can  greatly  improve  their 
accuracy.  Understanding  the  impact  of  these  processes  on  large  scale  circulations  can 
also  help  to  provide  better  ice  melt  forecasts,  as  well  as  higher  resolution  global  and 
regional  oceanographic  forecast  models  across  all  spectrums  within  oceanography. 
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Figure  2.  Diagram  illustrating  the  geographic  distribution  of  diffusive-convection, 
which  occurs  in  approximately  15%  of  the  world’s  oceans.  Unlike  salt  fingering, 
this  phenomenon  occurs  when  cold  fresh  water  is  located  above  warm  and  salty. 


Figure  3.  Diagram  illustrating  geographic  locations  of  salt  fingers,  which  occurs  in 

approximately  30%  of  the  world’s  oceans. 

B.  BEYOND  OCEANOGRAPHY 

By  extending  our  knowledge  of  oceanographic  double-diffusion  and  changing 
background  parameters  we  can  get  a  “feel”  of  how  double  diffusive  convection  changes 
in  other  environments. 

External  planetary  systems,  such  as  the  atmospheric  makeup  of  planets  within  our 
solar  system,  are  represented  by  a  much  lower  Prandtl  number  (Pr)  than  in  our  oceans, 
(Pr  =  10'5  verses  Pr  =  7).  Estimates  show  that  for  Jupiter’s  atmosphere  (radiation  driven), 
Pr  «  10'4  and  for  the  metallic  interior  (conduction  driven),  Pr  «  10'2  (Bagenal  et  al.„ 
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2004).  Decreasing  Pr  from  the  oceanographic  value  of  7,  allows  us  to  build  a  better 
understanding  of  how  double  diffusive  convection  behaves  within  these  foreign  bodies. 

The  interior  temperature  and  Helium  within  our  own  sun  and  other  stars  is  also  a 
natural  area  where  double  diffusive  convection  occurs  as  interior  temperature  and  Helium 
diffuse  at  separate  rates.  Metallic  lingers  are  theorized  in  these  systems  to  consist  of 
metallic  “blobs”  that  descend  similar  to  the  oceanographic  warmer  denser  fluid,  with  heat 
diffusing  much  quicker  than  the  heavy  metal  composition  within  the  fluid  (Vauclair 
2003).  Another  manifestation  of  active  salt  lingering  in  nature,  the  high  Prandtl  number 
within  magma  flows,  has  close  similarities  to  oceanographic  salt  lingers  (Tait  1989). 

Sugar-salt  interactions  have  been  modeled  in  the  laboratory  environment  starting 
with  Jevons  in  1857  with  pure  water  and,  hydrochloric  acid  and  white  sugar  (Jevons 
1857).  The  sugar-salt  solution  is  different  to  prior  cases,  in  both  the  Prandtl  number  and 
the  diffusivity  ratio.  However,  we  include  this  case  so  as  to  achieve  a  more  complete 
understanding. 

While  many  regimes  are  currently  beyond  the  modeling  capabilities  of  even  for 
the  most  powerful  super  computers,  in  this  study  we  strive  to  gain  a  quantitative 
understanding  within  the  accessible  limits  of  Pr,  x,  and  Rp  values.  Many  regimes  we  wish 
to  explore  are  outside  this  reach  and  we  are  pushing  the  limits  of  what  is  possible  with 
today’s  computer  processing  speeds. 

C.  EXTANT  MODELS  OF  EQUILIBRIUM  TRANSPORT 

A  fundamental  question  in  double-diffusive  convection  theory  concerns  the 
equilibrium  transport  of  temperature,  salinity,  chemical  tracers,  and  momentum.  The 
quantification  of  double-diffusive  fluxes  and  their  dependencies  on  the  background 
temperature  and  salinity  gradients — so-called  flux-gradient  laws — is  an  essential  step  in 
linking  the  microstructure  dynamics  with  larger  scales  of  motion.  This  problem  has 
motivated  numerous  laboratory  (Lambert  and  Demenkow  1972;  Griffiths  and  Ruddick 
1980)  and  field  (Schmitt  et  al.„  1987;  Schmitt  et  al.„  2005)  experiments,  as  well  as 
theoretical  (Stern  1969;  Kunze  1987)  and  numerical  (Shen  1995;  Stern  et  al.„  2001; 
Kimura  and  Smyth  2007,  2011;  Traxler  et  al.„  2011;  Stellmach  et  al.„  2011)  models. 
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While  the  linear  stability  theory  of  double  diffusion  (Stem  1960;  Baines  and  Gill  1969)  is 
well  developed  and  fully  understood,  it  does  not  explain  the  saturation  of  primary  double- 
diffusive  instabilities — the  problem  of  equilibrium  double-diffusive  transport  is 
complicated  and  principally  nonlinear. 

Numerous  attempts  have  been  made  to  deduce  flux-gradient  laws  from  first 
principles.  The  first  and  perhaps  the  most  influential  idea  was  introduced  by  Stern 
(1969),  who  suggested  that  the  linear  growth  of  salt  fingers  is  arrested  when  the  Stern 
number 

£  —  ~  PFsi\m 

~  v(af  -  [3 Sz) 

reaches  0(1).  FTdim  and  Fsdim  here  are  the  dimensional  temperature  and  salinity  fluxes;  T, 
and  S_  are  the  vertical  gradients;  a,  [1  are  the  thennal  expansion  and  haline  contraction 
coefficients  respectively;  v  is  the  molecular  viscosity.  Stern’s  physical  argument  is 
compelling.  When  A  exceeds  unity,  the  unbounded  salt  finger  system  becomes  unstable 
with  respect  to  collective  instability,  a  term  used  to  describe  the  spontaneous  excitation  of 
gravity  waves  by  salt  fingers.  Stem  suggested  that  collective  instability  could  dismpt  salt 
fingers,  thereby  arresting  the  development  of  primary  instabilities. 

At  first  glance,  various  pieces  of  indirect  evidence  seemed  to  validate  Stem’s 
hypothesis.  Oceanographic  measurements  (Hebert  1988;  Innoue  et  al.„  2008)  tend  to 
produce  0(1)  values  of  A.  Kunze  (1987)  gave  an  alternative  argument  in  support  of  the 
Stern  number  constraint.  He  showed  that  A  =  1  is  equivalent  to  the  Richardson  number 
(Ri)  of  %  based  on  scales  of  individual  salt  fingers.  Kunze  suggested  that  the  well-known 
criterion  for  the  instability  of  horizontal,  inviscid  and  nondiffusive  parallel  flows 
(Richardson  1920;  Howard  1961;  Miles  1961),  i.e.  Ri  <  %,  can  be  extended  to  viscous, 
diffusive  and  vertically  oriented  fingers. 

However,  closer  inspection  has  revealed  certain  inconsistencies  in  the  Stem- 
Kunze  hypothesis.  For  instance,  in  many  laboratory  experiments,  salt  and  sugar  replace 
heat  and  salt  as  buoyancy  components  (Stern  and  Turner  1969;  Lambert  and  Demenkow 
1972;  Griffiths  and  Ruddick  1980;  Taylor  and  Veronis  1996;  Krishnamurti  2003). 
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Typical  values  of  Stern  number  in  this  case  are  extremely  low.  Lambert  and  Demenkow 

•3 

(1972)  report  Stem  numbers  as  low  as  A  =  2-10'  ,  casting  some  doubt  on  the  generality  of 
the  Stern  number  constraint.  On  the  theoretical  side,  a  critical  advance  was  made  by 
Holyer  (1984),  who  performed  a  formal  linear  stability  analysis  for  vertical  steady  salt 
lingers.  She  discovered  direct,  relatively  small-scale  —  comparable  to  the  salt  finger 
width  —  secondary  instabilities,  which  typically  grow  much  faster  than  the  collective 
instability  modes.  Unlike  collective  instabilities,  Holyer  modes  grow  regardless  of  the 
(finite)  amplitude  of  salt  fingers  and  thus  regardless  of  the  specific  values  of  A. 
Numerical  simulations  (Shen  1995;  Radko  and  Stern  1999)  confirmed  Holyer’s  ideas  by 
demonstrating  that  the  equilibration  of  salt  fingers  occurs  by  means  of  the  nonlinear 
interaction  of  primary  instabilities  with  Holyer  modes.  Furthermore,  the  eddy 
diffusivities  of  heat  and  salt  from  simulations  (Stern  et  al.„  2001)  and  available 
observations  (St.  Laurent  and  Schmitt  1999)  monotonically  decrease  with  the  density 
ratio,  whereas  the  opposite  trend  was  expected  based  on  the  Stern-Kunze  constraint.  The 
most  recent  and  well-resolved  simulations  of  the  heat-salt  system  (Traxler  et  al.„  2011) 
indicate  that  as  the  density  ratio  increases  from  1.2  to  10,  the  Stern  number  in  the 
equilibrium  flow  reduces  by  more  than  two  orders  of  magnitude,  from  A  =  76  t o  ^4  =  0.6. 

Concerns  with  regard  to  the  Stern  number  constraint  have  motivated  the 
development  of  alternative  models  for  the  equilibrium  transport.  While  the  general 
analytical  description  of  double-diffusive  transport  is  still  lacking,  several  cases  have 
already  been  successfully  treated.  Weakly  nonlinear  instability  theory  describes  salt 
finger  patterns,  dynamics,  and  transport  characteristics  for  the  parameters  near  the  point 
of  marginal  instability  (Proctor  and  Holyer  1986;  Radko  and  Stern  1999,  2000;  Stern  and 
Simeonov  2004;  Radko  2010).  Promising  attempts  to  formulate  physically  based 
parameterizations  of  salt  finger  transport  include  the  mixing  length  model  (Stem  and 
Simeonov  2005),  a  double-diffusive  version  of  the  upper  bound  theory  (Balmforth  et  al., 
2006),  a  second-order  closure  approach  (Canuto  et  al.,  2008),  and  various 
phenomenological  models  (Shen  1995;  Merryfield  and  Grinder  1999;  Radko  2008).  We 
propose  what  appears  to  be,  to  date,  the  most  general  algorithm  for  estimating  the 
equilibrium  salt  finger  transport.  The  model  is  based  on  the  properties  of  secondary  salt 
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finger  instabilities  and  can  be  applied  to  a  variety  of  parameter  regimes,  including  small, 
moderate,  and  large  Prandtl  numbers,  and  a  wide  range  of  diffusivity  ratios. 


D.  OUR  APPROACH 

Our  story  is  a  new  variation  on  an  old  theme — a  theme  that  appears,  explicitly  or 
implicitly,  in  almost  all  theories  of  equilibrium  double-diffusive  convection  (Stern  1969; 
Holyer  1984;  Kunze  1987;  Stern  and  Simeonov  2005;  Radko  2010).  We  revisit  the  idea 
of  a  competition  between  linear  mechanisms  involved  in  the  growth  of  salt  fingers  and 
the  disruptive  action  of  their  secondary  instabilities.  At  the  same  time,  we  strive  to  avoid 
some  debatable  assumptions  of  earlier  models,  such  as  the  link  between  secondary 
instabilities  and  Stern  number.  The  competition  between  primary  and  secondary 
instabilities  is  concisely  phrased  in  the  growth  rate  balance  assumption: 

(1) 

where  X\  is  the  typical  growth  rate  of  linear  salt  fingers,  X2  is  the  growth  rate  of  their 
secondary  instabilities,  and  C  is  a  dimensionless  order  one  quantity  that  can  be  calibrated 
on  the  basis  of  simulations  or  experiments.  The  primary  growth  rate  k|  is  detennined  by 
the  background  gradients  of  temperature  and  salinity  and  by  molecular  characteristics 
(diffusivities  and  viscosity).  The  secondary  instability  X2  is  also  affected  by  these 
quantities,  but,  additionally,  it  depends  very  strongly  on  the  amplitude  of  primary  salt 
fingers.  Thus,  for  any  given  background  parameters  and  a  value  of  C,  the  growth  rate 
balance  Equation  (1)  implicitly  determines  the  equilibrium  amplitude  of  salt  fingers. 

The  physical  reasoning  behind  Equation  (1)  is  straightforward.  As  initially  weak 
salt  fingers,  emerging  from  random  small  scale  perturbations,  grow  in  time,  they  start  to 
develop  secondary  (Holyer  1984)  instabilities.  Unlike  the  primary  ones,  the  growth  rate 
of  secondary  instabilities  monotonically  increases  with  the  amplitude  of  fingers.  At  first, 
the  growth  of  secondary  instabilities  is  too  slow  to  inflict  any  significant  damage  on 
growing  primary  modes — the  evolution  of  small  amplitude  fingers  is  adequately  captured 
by  the  linear  theory.  However,  at  some  point,  the  growth  rate  of  secondary  instabilities 
exceeds  the  primary  growth  rates.  As  a  result,  the  secondary  instabilities  start  to  gain  in 
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magnitude,  rapidly  reaching  the  level  of  primary  modes.  Their  interaction  nonlinearly 
suppresses  the  growth  of  salt  lingers.  At  this  stage,  the  system  reaches  statistical 
equilibrium. 

In  this  study,  condition  Equation  (1)  is  used  as  the  basis  for  a  simple  algorithm  to 
determine  the  equilibrium  amplitude  of  salt  lingers.  The  growth  of  primary  modes  (ki)  is 
evaluated  from  the  linear  instability  theory  (Stern  1960;  Baines  and  Gill  1969).  For  any 
given  amplitude  of  salt  fingers,  the  growth  rate  of  secondary  instabilities  (X2)  can  be 
computed  using  elements  of  the  Floquet  theory  as  in  Holyer  (1984).  The  amplitude  of 
fingers  is  iteratively  adjusted  until  A\  and  I-2  satisfies  the  growth  rate  balance  Equation  (1) 

.  The  model  predictions  are  then  conveniently  expressed  in  terms  of  the  equilibrium  eddy 
fluxes  of  heat  and  salt. 

The  use  of  DNS  has  been  essential  to  the  success  of  our  research.  DNS  allows  us 
to  (i)  calibrate  our  theoretical  model,  (ii)  validate  our  reference  base  on  the  growth  rate 
balance  Equation  (1),  and  (iii)  expand  the  sensitivity  of  our  solutions  to  the  background 
parameters.  From  crude  calculations  of  man-hours  spent  on  this  project  it  would  take 
91+  years  (over  800,000  hrs)  to  run  these  calculations  on  a  single  computer  processor 
continuously.  Many  of  these  experiments  are  run  on  multiple  parallel  processors, 
sometimes  128  working  simultaneously  and  still  taking  longer  than  a  week  to  complete. 
As  computing  power  increases,  this  processing  time  will  decrease  and  more  robust 
calculations  and  data  can  be  collected. 

This  thesis  is  organized  as  follows.  In  Chapter  II,  we  present  DNS,  focusing  our 
inquiry  on  the  equilibration  of  primary  instability.  A  theoretical  model  of  equilibration  is 
presented  in  Chapter  III.  We  compute  the  temperature  and  salinity  fluxes  as  a  function  of 
density  ratio  for  the  oceanographic  case  (Pr=7,  r=0.01)  and  successfully  test  (Chapter  IV) 
the  theoretical  prediction  by  DNS.  In  Chapter  V,  we  also  explore  a  broader  parameter 
range,  including  the  low  Prandtl  number  regime,  relevant  for  astrophysical  applications, 
and  high  Prandtl  numbers,  relevant  for  magmatic  melts.  Each  case  is  compared  with  the 
corresponding  DNS.  We  summarize  and  draw  conclusions  in  Chapters  VI  while 
discussing  future  possibilities  in  Chapter  VII. 


9 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


10 


II.  PRELIMINARY  CALCULATIONS 


The  simplest  and  most  direct  approach  to  the  analysis  of  equilibrium  of  salt 
fingering  involves  Direct  Numerical  Simulation  (DNS).  We  used  a  de-aliased  pseudo- 
spectral  method  described  and  first  used  for  two-dimensional  (Stern  and  Radko  1998)  and 
three  dimensional  salt-fingering  simulations  (Radko  and  Stern  1999).  In  order  to 
expedite  the  simulations  the  numerical  algorithm  was  parallelized  using  a  Message 
Passing  Interface  (MPI)  algorithm  in  both  two  and  three  dimensions.  Two  dimensional 
model  runs  were  performed  on  Anastasia,  a  high-performance  computer  cluster  at  Naval 
Postgraduate  School  (NPS).  Three  dimensional  model  runs  were  perfonned  on  the 
shared  high-performance  computer  cluster  Hamming,  also  at  NPS,  and  Einstein,  a  U.S. 
Navy  Cray  XT5  high-perfonnance  computer  cluster.  The  resources  utilized,  varied  from 
4  processors  on  Anastasia  up  to  128  simultaneous  processors  on  Einstein  for  the  most 
highly  resolved  cases  in  heat  salt  throughout  Rp  values  of  1.1,  1.3,  ...,  2.9  (Tables  2 
through  6).  We  resolve  a  total  of  ten  fingers;  5  warmer,  more  dense  and  5  cooler,  less 
dense.  However,  in  order  to  resolve  the  salt  finger  scale,  these  runs  required  anywhere 
from  3  days  to  4  weeks  to  complete.  With  increased  computing  power,  this  time 
limitation  to  process  data  will  decrease  in  future  years.  In  turn,  this  will  allow  future 
graduate  students  the  flexibility  to  run  more  data  sets  over  longer  periods  and  thereby 
address  even  more  challenging  problems  in  double-diffusive  convection. 


Following  Radko  and  Stern  (1999),  we  separate  the  temperature  and  salinity  fields 
into  the  basic  state  (t>Sj ,  representing  a  uniform  vertical  gradient,  and  a  departure 

(T, S)  from  it.  The  two-dimensional  Boussinesq  equations  of  motion  are  expressed  in 


terms  of  T,S  and  nondimensionalized  using  /  = 


f  kT  v  V  k 


— ,  — ,  and  ^°V^T  as  the 


/  kT 


l 1 


\  §°~T-  j 

scales  of  length,  velocity,  time  and  pressure  respectively.  Here,  ( kr,ks )  denote  the 
molecular  diffusivities  of  heat  and  salt  and  po  is  the  reference  density  used  in  the 
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Boussinesq  approximation.  The  expansion/contraction  coefficients  are  incorporated  in 
(T,S)  and  aT zl  is  used  as  the  scale  for  both  temperature  and  salinity  perturbations, 
resulting  in 
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V- v  =  0, 


(2) 


where  v  =  (w,v,w)  is  the  velocity  vector.  This  system  is  unstable  with  respect  to  salt 
lingering  (Stem  1960)  for 

\<rp<\/t.  (3) 

The  key  nondimensional  numbers  governing  the  evolution  of  system  Equation  (2) 
are  the  Prandtl  number  pr  =  —  ,  the  diffusivity  ratio  r  =  — ,  and  the  background  density 

kT  kj. 

aT 

ratio  R  =  — =?-.  We  also  assume  that,  in  the  absence  of  large-scale  structures,  fluxes  are 

p  psz 

independent  of  the  nondimensional  parameters  related  to  the  domain  size  (e.g.,  the 
Rayleigh  number).  The  local  flux-gradient  laws  are  commonly  used  to  parameterize  the 
effects  of  salt  fingering  on  the  oceanic  circulation  and  successful  attempts  have  been 
made  to  validate  the  concept  of  an  unbounded  T-S  gradient  in  numerical  simulations 
(Radko  and  Stem  1999). 

To  gain  a  preliminary  understanding  of  the  mechanics  of  equilibration  (Rp,Pr,x), 

Equation  (2)  was  solved  numerically.  We  assume  triply  periodic  boundary  conditions  for 

T,  S,  p,  and  v  and  integrate  the  governing  equations  using  a  de-aliased  pseudo-spectral 

method  described  in  Radko  and  Stern  (1999).  In  the  following  calculation,  we  use 

parameters  relevant  in  the  oceanographic  (heat/salt)  context:  z=0.01  and  Pr=7.  The 

overall  density  ratio  is  RP=T.9,  which  is  representative  of  the  main  Atlantic  thermocline. 

The  size  of  the  computational  domain  (Lx  =  38.2,  Ly  =  38.2,  Lz  =  76.3)  is  equivalent  to 
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(5x5x10)  linearly  fastest  growing  finger  wavelengths  ( d ).  The  flow  was  resolved  by  a 
uniform  mesh  with  (Nx  x  Ny  x  Nz)  =  (256  x  256  x  512)  elements,  and  the  model  was 
initialized  from  rest  by  a  small-amplitude  random  computer-generated  initial  ( T,S ) 
distribution.  After  a  few  characteristic  growth  periods,  active  statistically  steady  double- 
diffusive  convection  was  established.  Figure  4a  shows  a  typical  instantaneous  (t=  20) 
temperature  field  in  the  initial  stage  of  linear  growth. 


Figure  4.  Equilibration  of  salt  fingers  in  the  numerical  experiment  with  Pr=7,  x  =  0.01, 
Rp=  1.9.  Three-dimensional  instantaneous  temperature  fields  are  shown  for  a) 
the  early  stage  of  linear  growth  at  t=20  and  b)  the  fully  equilibrated  state  at  t=  50. 
Red/green  color  corresponds  to  high  values  of  T  and  low  values  are  shown  in  blue 
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As  expected  from  the  linear  stability  theory,  the  most  rapidly  growing 
perturbations  take  the  form  of  the  vertically  uniform  salt  fingers  known  as  elevator 
modes.  Figure  4b  presents  the  fully  developed  equilibrium  regime  (t=50).  While  the 
previously  dominant  elevator  mode  is  still  visible,  it  is  now  comparable  in  magnitude  to 
the  secondary  instabilities,  which  act  to  twist  the  vertical  fingers,  adversely  affecting  their 
growth.  This  stage  is  characterized  by  the  irregular  transient  patterns  in  the  temperature 
field — a  dramatic  consequence  of  the  nonlinear  interaction  between  the  elevator  modes 
and  their  secondary  instabilities. 


Figure  5.  Time  record  of  the  temperature  (solid  curve)  and  salinity  (dashed  curve) 
fluxes  for  the  calculation  in  Fig  4.  Equilibrium  point  is  marked  by  red  line. 


The  transition  to  the  statistically  steady  regime  is  illustrated  in  Figure  5,  which 
shows  the  evolution  of  the  spatially  averaged  downward  temperature  and  salinity  fluxes 
(Ft,Fs)  =  (-wT,-wS)  -  the  sign  convention  is  dictated  by  considerations  of  convenience. 
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Note  that  the  heat  flux  in  our  system  of  nondimensionalization  is  FT=Nu- 1 ,  where  Nu 
is  the  Nusselt  number — the  ratio  of  the  total  heat  flux  to  the  (much  weaker)  molecular 
diffusive  flux.  The  initial  stage  of  the  simulation  (t<30)  is  characterized  by  the 
exponential  growth  of  fluxes,  which  is  followed  by  their  equilibration.  After 
equilibration,  the  intensity  of  salt  fingering  remains  statistically  steady. 

The  nondimensional  fluxes  ( FT ,  Fs )  can  be  converted  to  the  dimensional  eddy 
diffusivities  of  heat  and  salt  (KTdim,K Sdim)  using 

^ t dim  =  FTkT,  KSiim  =  FskTRp.  (4) 

The  time  mean  T-S  fluxes,  averaged  over  the  second  half  of  the  experiment  in 
Figure  5,  are  FT  =40.88  ,  Fs  =75.88,  which  corresponds  to  the  dimensional  diffusivities 
of  heat  and  salt  of 

KTAim  =5.6-10-6mV1,  KSdim  =  2.0  •  10-5m  V1  (5) 

These  values  are  broadly  consistent  with  the  observational  estimates  (St.  Laurent 
and  Schmitt,  1999)  for  the  density  ratio  of  Rp=  1.9. 

In  Figure  6,  we  plot  typical  horizontal  (top  panel)  and  vertical  (lower  panels) 
sections  of  temperature  (left)  and  salinity  (right)  for  the  fully  developed  equilibrium 
regime.  Temperature  and  salinity  patterns  are  correlated  and  both  reveal  vertically 
elongated — but  far  from  uniform — horizontally  isotropic  salt  fingers.  Salinity  sections 
are  characterized  by  very  fine,  relative  to  temperature,  spatial  scales,  which  is  a  direct 
consequence  of  the  low  diffusivity  ratio. 
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Figure  6.  Horizontal  (upper  panels)  and  vertical  (lower  panels)  sections  of  temperature 
(left)  and  salinity  (right)  for  a  typical  fully  equilibrated  state  of  the  calculation  in 

Fig  5. 
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III.  THEORETICAL  MODEL 


Guided  by  physical  reasoning  (Section  D  of  Chapter  I)  and  by  direct  numerical 
simulations  (e.g.,  Chapter  II),  we  now  attempt  to  predict  the  equilibrium  level  of  salt 
fingers  from  the  growth  rate  balance  Equation  (1).  The  growth  rate  equation  for  the 
primary  salt  finger  instabilities  (Stem  1960)  in  our  nondimensional  units  reduces  to 

A3  +  [l  +  r  +  Pr]A :2A2  +  [(r  +  Pr+  rPr)/t4  +  Pr  (l  -  Rp'  )]  A  +  r  Pr  k6  -Pr  k2  (i?;1  -  r)  =  0,  ^ 

where  k  is  the  horizontal  wavenumber  of  the  vertically  uniform  elevator  modes.  For  our 
purpose  it  is  not  necessary  to  consider  inclined  salt  fingers  since  the  largest  growth  rates 
are  attained  by  the  vertically  oriented  elevator  modes  (e.g.,  Baines  and  Gill  1969).  For 
each  (Rp,Pr,r),  we  compute  the  largest  growth  rate  (A|)  by  maximizing  the  solutions  of 
Equation  (6)  with  respect  to  k  as  in  Schmitt  (1979,  1983).  Thus,  A|  is  uniquely 
determined  by  the  governing  parameters: 

4=4(*„,Pr,r).  (7) 

While  the  maximal  growth  rates  of  primary  instabilities  are  identical  in  two  (x,z) 
and  three  (x,y,z)  dimensions,  the  secondary  instabilities  differ  substantially.  Therefore, 
the  algorithms  for  computing  An  in  two  and  three  dimensions  will  be  discussed  separately. 


A.  TWO-DIMENSIONAL  FORMULATION 

In  two  dimensions,  the  governing  equations  are  conveniently  written  in  terms  of 

the  streamfunction  1//  ,  such  that  (u,w)  =  (  -- — j ,  resulting  in 

V  dz  dx  J 


dT  dw  , 

—  +  T)  H — —  =  V2T 

dt  dx 

dS  \  dw  2 

—  +  J{y/,S)  + - -  =  zV5 

dt  R  dx 


(8) 


d  ■,  . 

—  VV  +  J{y/,V~y/)  =  Pr 
dt 


A(r_5)  +  vy 

dx 


Our  next  step  is  to  use  Equation  (8)  to  examine  the  stability  of  the  basic  state 
representing  the  elevator  modes: 
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T  =  At  sin (kx) 

<  S  -  As  sin(foc)  (9) 

\j7  =  Aif/  cos  (kx) 

It  should  be  mentioned  that  the  secondary  stability  problem  can  be  formulated  in 
different  ways.  Holyer  (1984)  has  chosen  the  wavenumber  k  in  Equation  (9)  to 
correspond  to  the  solution  with  zero  primary  growth  rate  (k0).  This  restriction  defines  a 
well-posed  stability  problem  with  the  regular  basic  steady  state.  There  is  only  one  reason 
to  question  this  formulation.  It  is  apparent  from  DNS  (Stern  et  al.„  2001;  Traxler  et  al.„ 
2011)  and  even  indicated  by  some  oceanographic  field  measurements  (Gargett  and 
Schmitt,  1982)  that  the  horizontal  wavenumber  of  dominant  fully  developed  salt  lingers 
is  more  accurately  approximated  by  kmax,  the  wavenumber  corresponding  to  the 
maximum  growth  rate,  than  by  k0.  Use  of  k=kmax  in  Equation  (9);  however,  results  in  the 
time-dependent  basic  state,  and  therefore  the  conventional  methods  of  stability  analysis 
cannot  be  applied  directly.  This  problem  is  not  uncommon  and  not  insurmountable. 
Notable  examples  of  analyzing  such  flows  come  from  studies  based  on  the  Kolmogorov 
model — the  sinusoidal  parallel  flow  in  viscous  fluid  (Sivashinsky  1985;  Manfroy  and 
Young,  1999  2002;  Balmforth  and  Young  2002,  2005).  Since  the  unforced  viscous 
parallel  flow  would  inevitably  decay,  the  Kolmogorov  model  circumvents  this  problem 
by  introducing  artificial  forcing  in  the  momentum  equation  that  maintains  the  steady 
state.  Variations  on  the  same  principle,  often  referred  to  as  the  quasi  steady  state 
approximation  or  a  “frozen  flow”  method,  have  been  successfully  applied  to  numerous 
stability  problems  (e.g.,  Lick  1964;  Robinson  1976)  including  salt  fingers  (Kimura  and 
Smyth  20 11). 

A  choice  has  to  be  made  at  this  point  between  the  approach  of  a  purist,  insisting 
on  formal  mathematical  consistency,  and  that  of  a  practically-oriented  researcher,  more 
concerned  by  the  consistency  with  patterns  realized  in  nature.  We  have  examined  both 
models,  k=ko  and  k=kmax,  and  the  results  are  qualitatively  similar.  The  prediction  of  the 
k=kmax  model  is  closer  to  the  corresponding  DNS,  and  therefore  we  limit  the  following 
discussion  to  the  configuration  based  on  the  fastest  growing  solution.  We  assume  that  in 
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the  regime  close  to  equilibrium  the  amplitude  Equation  (9)  is  maintained  at  a  quasi¬ 
steady  level  for  k=kmax  and  separate  the  solution  into  the  dominant  basic  state  and  small 
perturbation: 

{T,S,yf)  =  {f,S,ijf)+{T',S',yf').  (10) 


Equation  (8)  is  linearized  as  follows: 


dT'  ,  ,  .  .ar  ,  ,  /(  .dur'  dy/’  „2rnl 

dt  V  dz  T  dz  dx 

dS'  dS'  .  ,  ,,  .  du A  1  dy/'  _2o, 

- A  k  sin  (kx) - A,k  cos  (kx) - h - =  rV~S 

dt  v  dz  s  dz  R  dx 


(ID 


(5  2  ,  .  ,  .  .dV  t//'  3  dy/' 

—  Vi//  - Av/ksm(kx)  - 4/  sin(^r)^  =  Pr 


—(T'-S')  +  vY 
dx 


To  examine  the  stability  properties  of  the  linear  system  Equation  (11),  we  use  the 
Floquet  technique  in  which  the  perturbation  is  sought  in  the  following  form: 


r J 

N 

(  t\ 

1  n 

S' 

=  exp {ifkx  +  imz  +  At)  ^ 

K 

y ) 

n=-N 

exp(inkx) , 


(12) 


where  k  is  the  growth  rate,  m  is  the  vertical  wavenumber,  and/is  the  Floquet  coefficient, 
which  controls  the  fundamental  wavelength  in  x.  Substituting  Equation  (12)  into  the 
linear  system  Equation  (11)  and  collecting  the  individual  Fourier  components  allows  us 
to  express  the  governing  equations  in  the  matrix  form: 


ii 

'Axj 

(13) 

where 

h  —  (T_N,S_N,y/_N,T_N+l,S_N+x,ii/_N+x,  ...  ,Tn,Sn,i//n)  , 

(14) 

and  A 

is 

the  square  matrix  whose  elements  are 

functions 

of[k,m,f,Rp,Vx,r,N,AT,As,Av). 

The  growth  rates  of  the  normal  modes  of  the  linear  system  Equation  (11) 
correspond  to  the  eigenvalues  of  A  .  For  each  set  of  governing  parameters,  we  determine 
the  eigenvalue  with  the  maximum  real  part,  which  represents  the  fastest  growing  mode. 
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The  value  of  this  growth  rate  is  then  maximized  with  respect  to  m  and  f  Note  that  the 
symmetry  and  periodicity  properties  of  our  system  are  such  that  the  Floquet  coefficient/ 
needs  to  be  varied  only  within  the  range  0  <  f  <  0.5 . 


The  coefficients  (AT^4S,AV)  of  the  normal  mode  Equation  (9)  are  connected,  for 
k=kmax,  by  the  following  relations: 


AT(A  +  k2  )  =  k  A 

T  \  1  max  J  max  y/ 


(15) 


and  therefore  the  amplitude  vector  (AT,AS^4V)  can  be  immediately  determined  from  any 
component — for  instance  from  the  temperature  amplitude  (At).  As  described  below,  the 
solution  rapidly  converges  with  increasing  resolution  (N—> oo).  Thus,  the  fastest  growing 
secondary  instability  is  determined  by  four  parameters: 

A7=A7(Rp,Pr,r,AT).  (16) 


Our  next  step  is  to  solve  the  growth  rate  balance  Equation  (1).  For  each  set  of 
governing  parameters  (Rp,Pr,x),  we  determine  the  primary  growth  rate  /  using  Equation 
(6).  Suppose  now  that  the  value  of  C  in  Equation  (1)  is  known.  Then,  by  varying  the 
amplitude  of  the  normal  mode  Aj  we  can  readjust  Equation  (16)  until  the  growth  rate 
balance  is  satisfied.  The  dependence  AMt)  takes  the  form  of  a  monotonically  increasing, 
nearly  linear  relation.  For  Aj  =  0,  the  secondary  instability  problem  becomes  identical  to 
that  of  primary  instability  and  therefore  X2  =  Xt  ■  A  corollary  of  this  observation  is  that 
solutions  of  Equation  (1)  exist  only  for  Ol,  as  indicated  by  the  schematic  in  Figure  7. 
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Figure  7.  Illustration  of  the  algorithm  for  predicting  the  equilibrium  amplitude  of  salt 
lingers.  The  growth  rate  of  the  secondary  instability  (X2)  monotonically  increases 
with  the  amplitude  of  salt  fingers  ( Ar ).  We  hypothesize  that  the  equilibrium  (A*) 

is  reached  when  X2=C  Xi. 


An  iterative  procedure  for  solving  the  growth  rate  balance  Equation  (1)  was  coded 
in  MAPLE  and  it  typically  required  7-8  iterations  for  the  model  to  converge,  within  a 
negligible  relative  error  of  10'5,  to  the  sought-after  solution^  =  T*  .  The  vertical  fluxes 
are  then  reconstructed  using  Equation  (9)  and  Equation  (15)  as  follows: 


2  f? 


(4  +  ^max  ) 


(17) 
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Figure  8.  Heat  flux  in  the  theoretical  two-dimensional  model  for  various  values  of  C. 

In  Figure  8,  the  predicted  heat  and  salt  (Pr=7,  x=0.01)  fluxes  are  plotted  as  a 
function  of  density  ratio  (Rp)  for  various  values  of  C.  Fluxes  increase  with  C  and  rapidly 
reduce  with  Rp.  The  calculations  in  Figure  8  were  perfonned  for  N= 32.  However,  the 
specific  number  of  Fourier  harmonics  involved  in  the  Floquet  calculation  is  not 
particularly  important.  The  rapid  convergence  of  the  solution  with  increasing  N  is 
indicated  in  Table  1.  The  very  crude  truncation  with  N=  2  results  in  19%  error,  whereas 
the  N=4  error  is  less  than  2%. 
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N 

Ft 

Fs 

Ft  —Ft 

T  T  N=32 

err  = - . — 

FT\ 

T  W=32 

2 

29.55387 

50.23965 

0.1899 

4 

25.28658 

42.98554 

0.0181 

8 

24.83962 

42.22573 

0.0001 

16 

24.83668 

42.22074 

err  <  1CT5 

32 

24.83668 

42.22074 

NA 

Table  1 .  Convergence  of  the  Floquet  based  algorithm  in  a  two-dimensional 
calculation  for  (Rp,Pr,x)  =  (1.9,7,0.01).  As  the  number  of  spectral  harmonics  ( N ) 
increases,  the  heat  flux  (Ft)  rapidly  approaches  its  limiting  (N— »oo)  value. 


B.  THREE-DIMENSIONAL  FORMULATION 

Numerical  simulations  (e.g.,  Stern  et  al.„  2001)  reveal  significant  differences,  in 
magnitude  and  dynamics,  between  two-  and  three-dimensional  salt  fingers.  Three- 
dimensional  fluxes  are  systematically  higher,  by  a  factor  of  2-3,  than  their  2D 
counterparts.  In  this  section,  we  attempt  to  extend  the  growth  rate  balance  theory  to  three 
dimensions.  In  addition  to  the  natural  interest  in  three-dimensional  dynamics — after  all, 
we  live  in  a  three-dimensional  world — we  are  also  motivated  by  the  expectation  that  the 
results  will  help  to  rationalize  the  difference  in  the  intensity  of  salt  lingers  in  2D  and  3D. 
Since  the  following  theory  is  analogous  to  the  two-dimensional  case  (Section  A  of 
Chapter  III),  it  is  presented  in  abbreviated  form. 

The  essential  difference  between  two-  and  three-dimensional  dynamics  stems 
from  the  more  complicated  structure  of  the  salt  linger  elevator  mode  in  3D: 
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T  =Arex  p(Att)0(x,y) 

<  S  -  As  Qxp(\t)(j)(x,  y)  (18) 

w  =  Aw  exp (\t)(j>(x,y), 


where  the  horizontal  planform  function  <f>  satisfies  the  Helmholtz  equation  V2H</>  +  k 2<j>  -  0  . 
In  3D,  identical  growth  rates  can  be  attained  by  various  horizontal  planforms.  For 
instance,  even  if  we  limit  our  analysis  to  rectangular  planforms  cf>  =  cos(kxx)cos(kvy) ,  all 


possible  configurations  with 


+  K  =  k, 


2 

max 


(19) 


yield  the  maximum  growth  rates  (A|)  and  therefore  are  a  priori  equivalent.  The 
degeneracy  of  linear  theory  with  respect  to  the  planform  selection  is  well  known  and  can 
be  resolved  by  invoking  fundamentally  nonlinear  considerations.  Proctor  and  Holyer 
(1986)  and  Radko  and  Stern  (2000)  examined  the  problem  for  a  “bounded”  configuration, 
in  which  finger  layer  was  limited  by  rigid  horizontal  surfaces — the  setup  analogous  to  the 
classical  Rayleigh  convection  problem.  The  former  (latter)  study  suggested  the  tendency 
of  the  planform  to  evolve  spontaneously  to  two-dimensional  rolls  (square  cells).  For  the 
unbounded  regime,  some  guidance  regarding  the  planform  selection  was  provided  by 
diagnostics  of  numerical  simulations  in  Radko  and  Stern  (1999),  which  indicate  the 
preferred  planform  is  neither  square  nor  a  roll,  but,  rather,  some  combination  thereof, 
perhaps  best  described  by  the  planfonn  function  of  the  following  type: 

(/>  =  cos(A:maxx)  +  ^-cos(&max_v) ,  (20) 

which  was  used  in  all  calculations  described  below.  It  should  be  mentioned,  however, 
that  for  our  purpose — determination  of  the  equilibrium  transport  characteristics — the 
planform  choice  turns  out  not  to  be  particularly  significant,  with  the  square  cell  model 

performing  only  slightly  worse  than  Equation  (20).  1 


1  We  also  note  in  passing  that  sheared  environments  favor  formation  of  salt  sheets  aligned  in  the 
direction  of  the  background  current  (Linden,  1974;  Kimura  and  Smyth,  2007),  in  which  case  salt  finger 
dynamics  become  similar  to  the  2D  case  discussed  in  Sec.  3a.  It  is  perhaps  ironic  that  since  large-scale 
shears  are  ubiquitous  in  the  ocean,  salt  fingers  may  be  better  represented  by  two-dimensional  than  by  three- 
dimensional  simulations. 
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The  calculation  of  secondary  growth  rates  for  the  square  salt  lingers  Equation  (20) 
proceeds,  as  previously,  by  the  method  based  on  the  Floquet  theory.  The  governing 
equations,  Equation  (2),  are  linearized  with  respect  to  the  basic  state  Equation  (18); 
which  is  assumed  to  be  maintained  in  a  quasi-steady  state,  and  the  solution  is  sought  in 
the  following  form: 


nx,ny 

S' 

S*„ny 

u 

.J 

N  N 

Un 

=  exp (ifxkx  +  ifyky  +  imz  +  At)  ^  ^ 

nx,ny 

V 

nx=-Nny=-N 

V 

nx,ny 

W 

W 

nx  ,ny 

exp  (inrkx  +  invky ) , 


(21) 


where  (fx,fy)  are  the  Floquet  factors  in  x  and  y.  Substituting  Equation  (21)  in  the 
linearized  governing  equations  and  collecting  the  individual  Fourier  components  reduces 
the  stability  problem  to  matrix  form  Equation  (13).  Maximizing  the  growth  rates  with 
respect  to  (fx,fy,m),  we  predict  the  largest  growth  rates  of  secondary  instabilities  (m)  as  a 
function  of  (Rp,Pr,x,y4T)-  Finally,  for  any  given  C,  we  search  for  the  amplitude  Ar  =  A* 
that  satisfies  the  growth  rate  balance  Equation  (1).  The  vertical  fluxes  are  then 
reconstructed  as  follows: 

FT=^T=5-4(\+k2m  „) 

'  F  -o  5  4(4+*L)2  (22) 

=  wS  = - - - — — 

[  (4+r*L) 


and  plotted  in  Figure  9  as  a  function  of  density  ratio  (Rp)  for  various  values  of  C. 
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p 

Figure  9.  Fleat  flux  in  the  theoretical  three-dimensional  model  for  various  values  of  C. 

While  the  patterns  of  3D  (Figure  9)  and  2D  (Figure  8)  fluxes  are  qualitatively 
similar — in  both  cases  fluxes  rapidly  reduce  with  Rp  and  increase  with  C — we  note  that 
the  3D  fluxes  are  consistently  higher.  The  reason  is  attributed  to  the  fact  that  the  2D 
basic  state  Equation  (9)  is  more  unstable  than  its  3D  counterparts  Equation  (18)  and 
Equation  (20).  Thus,  the  amplitude  of  salt  fingers  has  to  be  significantly  higher  in  3D  to 
produce  secondary  instabilities  of  the  same  strength  as  in  2D,  which  rationalizes  the 
behavior  noticed  in  the  numerical  simulations  of  Stern  et  ah,  (2001). 
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IV.  OCEANOGRAPHIC  (HEAT-SALT)  CASE:  DNS  ANALYSIS 


In  this  chapter,  we  present  numerical  simulations  of  oceanographic  double- 
diffusive  convection,  warm  saltier  water  over  colder  fresher  water,  with  density  ratio  (Rp) 
values  1.1,  1.3,  ...,  2.9  (Table  2).  We  have  initiated  these  experiments  by  computer¬ 
generated  small  amplitude  temperature  and  salinity  distribution.  Although  double- 
diffusive  interactions  may  lead  to  larger  scale  phenomena,  such  as  thermohaline 
staircases,  thermohaline  intrusions,  and  internal  waves  (Traxler  et  ah,  2011),  this 
experiment  only  focuses  on  10  individual  fingers  (5  warm/salty  and  5  cool/fresh 
initially).  The  vertical  salinity  (Fs)  and  temperature  (Ft)  fluxes  were  recorded  as  a 
function  of  time  and  their  statistical  averages  were  analyzed  as  a  function  of  density  ratio. 
Simulations  of  this  nature  have  only  recently  become  possible,  as  computer  technology 
has  exceeded  the  petaflop  barrier.  This  has  allowed  us  to  explore  other  parameter  regimes 
besides  just  heat  and  salt  interactions  as  discussed  in  Chapter  V. 

A.  CALIBRATION  OF  THE  THEORETICAL  MODEL 

A  limitation  of  the  theoretical  model  in  Chapter  III  is  related  to  its  inability  to 
determine  the  constant  C  internally.  Our  theory  assumes  that  C  is  comparable  to  unity 
but  does  not  provide  an  exact  value.  In  order  to  calibrate  C,  we  now  turn  to  DNS.  For 
the  heat/salt  case  (Pr=7,T=0.01)  discussed  in  this  section,  the  analysis  is  based  on  10  two- 
dimensional  and  10  three-dimensional  simulations  conducted  for  Rp  =  1. 1,1.3,.  ..,2.9. 

Figure  10  presents  fluxes  from  a  series  of  two-dimensional  calculations,  plotted 
along  with  the  theoretical  prediction  from  the  foregoing  growth  rate  balance  theory 
(Chapter  II)  with  C=4.3.  In  all  simulations  we  used  (NX,NZ)=(256,5 1 2)  grid  points  and 
the  size  of  the  computational  domain  (LX,LZ)=(38.2,76.3)  equivalent  to  5  x  10  linearly 
fastest  growing  finger  wavelengths  ( d ).  The  variation  of  heat  (Figure  10a)  and  salt 
(Figure  10b)  fluxes  with  density  ratio  follows  the  theoretical  pattern  very  closely.  The 
ability  of  our  model,  containing  only  one  adjustable  parameter,  to  capture  the  details  of 
the  numerical  simulations  to  such  an  extent  lends  credence  to  the  proposed  condition 

Equation  (1)  as  the  most  relevant  physical  factor  constraining  growth  of  salt  fingers. 
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Figure  10.  Comparison  of  the  theoretical  two-dimensional  model  for  C=4.3  (solid  curve) 
with  DNS  (plus  signs).  Heat  and  salt  fluxes  are  shown  in  a)  and  b )  respectively. 

Figure  1 1  presents  the  corresponding  three-dimensional  calculations.  The 
resolution  of  three-dimensional  DNS  matched  the  two-dimensional  runs: 
(Nx,Ny,Nz)=(256,256,512)  mesh  was  used  to  resolve  5x5x10  fastest  growing  finger 
wavelengths.  The  theoretical  calculation  for  C=2.7  is  superimposed  on  the  numerical 
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data,  revealing  their  mutual  consistency — our  theory  predicts  the  right  order  of  magnitude 
and  the  pattern  of  fluxes  as  a  function  of  Rp.  The  relevant  value  of  the  adjustable 
coefficient  C=2.7  in  3D  is  comparable  to — but  differs  from — the  2D  value  of  C=4.3.  The 
(limited)  variation  in  C  can  be  readily  attributed  to  distinct  equilibrium  dynamics  of 
three-dimensional  and  two-dimensional  salt  fingers  in  the  unbounded  model  (Radko  and 
Stern  1999). 


Figure  1 1 .  The  same  as  in  Figure  10  but  for  the  three-dimensional  model  with  0=2.1 . 
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B. 


OBSERVED  PATTERNS 


T(x,y)  S(x,y) 


Figure  12.  Horizontal  distribution  of  temperature  and  salinity  for  fully  equilibrated  state. 

Oceanographic  values  of  Pr  =  7,  x  =  0.01,  Rp  =  1.9. 

The  horizontal  and  vertical  structure  of  individual  fingers  can  be  clearly  seen  from 
Figures  12  and  13,  respectfully.  It  is  interesting  to  note  that  the  temperature  field  is 
characterized  by  much  larger  scales  than  that  of  salinity.  It  is  also  evident  that  while  both 
the  temperature  and  salinity  distributions  are  spatially  correlated  there  is  more  variability 
within  the  salinity  field.  This  latter  feature  is  due  to  the  fact  that  salinity  diffuses  on 
much  smaller  scales  (100  times  slower)  than  temperature. 
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Figure  13.  Vertical  distribution  of  temperature  and  salinity  for  the  fully  equilibrated  state. 

Oceanographic  values  of  Pr  =  7,  x  =  0.01,  Rp  =  1.9. 
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c. 


MIXING  CHARACTERISTICS 


Equilibrium  Transport  -  Heat/Salt  -  Rhol  .9 


Figure  14.  Temperature  and  salinity  fluxes  as  a  function  of  time  for  the  heat-salt  case.  Ft 
and  Fs  values  are  calculated  after  equilibrium  is  reached. 

Temperature  and  salinity  fluxes  for  a  typical  North  Atlantic  oceanographic  case 
(Pr  =  7,  x  =  0.01,  Rp=  1.9)  are  plotted  in  Figure  14.  Note  that  fluxes  in  Figure  14  are 
given  in  dimensional  units  and  the  sign  convention  reflects  downward  transport  of  heat 
and  salt  (unlike  Figure  5).  Our  model  starts  using  random  noise  until  equilibrium  is 
reached  and  the  data  for  oceanographic  applications  are  taken  from  the  fully  equilibrated 
stage  (lower  panel  of  Figure  14). 

Using  data  similar  to  that  in  Figure  14  for  a  range  of  Rp,  we  can  calculate  a  mean 
of  our  fluxes  after  equilibrium  is  reached.  This  corresponds  to  typical  conditions  in  the 
Atlantic  thermocline  for  a  range  of  density  ratio  values  (Rp).  These  values  are  included  in 
Table  2,  along  with  the  dimensional  temperature  and  salinity  diffusivities  (Kjdim  and 
Ksdim)  which  are  of  the  same  order  of  magnitude  expected  in  the  ocean.  Also  included  in 

F 

Table  2  are  Buoyancy  Flux  (Fs  -  Ft)  and  Flux  Ratio  (y  =  — )  values  for  the  heat  salt  case. 
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Heat/Salt 

256x256x512 

Rp 

ft 

Fs 

KTdim(mA2/s) 

Ksdim(mA2/s) 

Buoyancy 

Flux 

Ratio 

l.i 

-253.3917 

-369.8733 

-3.471E-05 

-5.574E-05 

116.4816 

0.685077025 

1.3 

-114.1179 

-186.4594 

-1.563E-05 

-3.321E-05 

72.3415 

0.61202546 

1.5 

-70.2319 

-121.8825 

-9.622E-06 

-2.505E-05 

51.6506 

0.576226284 

1.7 

-49.5589 

-90.2265 

-6.790E-06 

-2.101E-05 

40.6676 

0.54927211 

1.9 

-40.8760 

-75.8794 

-5.600E-06 

-1.975E-05 

35.0034 

0.538696932 

2.1 

-35.0073 

-65.8659 

-4.796E-06 

-1.895E-05 

30.8586 

0.531493535 

2.3 

-28.9132 

-56.0535 

-3.961E-06 

-1.766E-05 

27.1403 

0.515814356 

2.5 

-27.3972 

-52.2506 

-3.753E-06 

-1.790E-05 

24.8534 

0.524342304 

2.7 

-23.0670 

-45.2340 

-3.160E-06 

-1.673E-05 

22.1670 

0.509948269 

2.9 

-22.7098 

-44.1900 

-3.11  IE-06 

-1.756E-05 

21.4802 

0.51391265 

Table  2.  Heat  -  Salt  temperatures  (Ft)  and  salinity  (Fs)  fluxes  for  a  range  of  density 
ratios,  Rp.  Also  shown  are  the  dimensionalized  temperature  and  salinity 
diffusivities  (KTdim,  Ksdim),  Buoyancy  Flux  (Fs-FT),  and  Flux  Ratio  (FT/Fs). 


While  development  of  explicit  parameterizations  of  double  diffusion  for  large- 
scale  ocean  modeling  is  outside  of  our  immediate  goal,  it  could  prove  beneficial  for  more 
applied  studies  to  note  that  the  numerical  data  in  Figure  1 1  can  be  accurately 
approximated  by  the  following  algebraic  expressions: 


Fs  =  °s  +bs,  (as,bs)  =  (135.7,-62.75) 

r  =  aY  exp (~brRp)  +  cy,  (ay,by,cy)  =  (2.709,2.513,0.5128) 
Ft  -  yFs, 


(23) 


which  are  indicated  in  Figure  15  by  the  red  curve.  The  structure  of  the  expression  for 
Fs(Rp)  is  suggested  by  theoretical  arguments  in  Radko  (2008),  whereas  the  expression  for 
the  flux  ratio  y(Rp)  represents  a  purely  empirical  exponential  fit  to  the  numerical  data. 
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Figure  15.  The  best  fit  for  Heat  Salt  (Pr=7;  x=0.01)  utilizing  cftool  in  MATLAB  to  solve 

Equation  (23). 
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V.  GENERALIZATIONS 


In  this  section,  we  go  beyond  the  oceanographically  motivated  heat-salt 
(Pr,x)=(7,0.01)  case  and  test  the  generality  of  our  conclusions  by  exploring  a  broader 
parameter  range.  The  discussion  in  this  section  is  based  on  the  comparison  of  four 
representative  series  of  three-dimensional  numerical  experiments  with  corresponding 
theoretical  results: 

•  Case  1 :  Low  Prandtl  number  (Pr,x)=(0. 1 ,0.0 1 ); 

•  Case  2:  High  Prandtl  number  (Pr,x)=(100, 0.01); 

•  Case  3:  High  diffusivity  ratio  (Pr,x)=(7, 0.1); 

•  Case  4:  High  Prandtl  number  and  High  diffusivity  ratio(Pr,x)=(1000, 1/3); 

Unfortunately,  DNS  for  more  extreme  parameter  values  (Pr— »co,  Pr— >0,  x— >0) 
become  computationally  prohibitive  in  3D.  However,  our  selection  represents  a 
substantial  fraction  of  the  numerically  accessible  range.  In  all  cases,  we  performed  a 
series  of  calculations  with  Rp  =  1.1, 1.3,..., 2. 9  in  computational  domains  corresponding  to 
5  x  5  x  10  fastest  growing  finger  wavelengths,  recording  the  average  equilibrium  fluxes. 

Comparing  Case  1  with  our  baseline  heat-salt  configuration  allows  us  to  examine 
specifics  of  the  low  Prandtl  number  regime,  which  is  relevant  for  astrophysical 
applications.  Likewise,  differences  between  Case  2  and  heat-salt  provide  some  glimpse 
into  the  large  Prandtl  number  regime,  exemplified  by  double-diffusion  in  magma 
chambers.  Case  3  differs  from  the  heat-salt  calculation  by  diffusivity  ratio,  which  also 
varies  dramatically  between  various  applications — from  x  =  0.8  for  heat/humidity  in  the 
atmosphere  to  x~T0'6  in  stellar  interiors.  Finally,  Case  4  (Sugar  Salt  parameters)  is  a 
combination  of  an  increase  in  Prandtl  number  and  diffusivity  ratio  and  behaves 
differently  than  the  majority  of  runs  throughout  our  experiments,  but  is  referenced  to 
show  how  these  differences  in  parameters  affect  results. 
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Density  Ratio 

Rp  =  1.9 

Pr 

X 

Ft 

Fs 

Heat- Salt 

7 

0.01 

-40.8760 

-75.8794 

Case  1 

0.1 

0.01 

-1.9140 

-4.6501 

Case  2 

100 

0.01 

-220.1848 

-364.8436 

Case  3 

7 

0.1 

-39.3424 

-60.2595 

Case  4 

1000 

1/3 

-769.4803 

-885.3976 

Table  3.  Dependence  of  the  equilibrium  T-S  fluxes  on  the  molecular  characteristics 
(Pr,x)  in  three-dimensional  simulations  with  Rp=1.9. 


A.  OBSERVED  PATTERNS 

We  can  gain  a  qualitative  understanding  of  the  effects  of  parameter  changes  by 
reviewing  Figures  16  to  19.  We  note  that  for  each  of  these  figures,  the  equilibrium  point 
has  already  been  reached. 
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Case  1 


Case  2 


Pr=7;t=0.1  Pi=1000;  x=l/3 

Figure  16.  DNS  model  of  horizontal  temperature  (highest  rate  of  diffusion)  in  (x,  y) 

coordinates  for  Cases  1  through  4. 

Taking  a  horizontal  slice  of  the  temperature  field  (Figure  16)  allows  us  to  further 
explore  the  similarities  and  changes  in  fingers  caused  by  differences  in  background 
parameters.  It  can  be  seen  that  there  exists  decreasing  dependence  between  Prandtl 
number  and  horizontal  finger  scales. 
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Case  1 


Case  2 


Pr=7;t=0.1 


I 


I 


Figure  17.  DNS  model  of  horizontal  salinity  (least  rate  of  diffusion)  in  (x,y)  coordinates 

for  Cases  1  through  4. 


At  first  glance,  a  review  of  the  horizontal  sections  of  salinity  fields  (Figure  17) 
does  not  show  as  clear  cut  of  a  correlation  as  that  of  temperature  fields.  By  breaking  the 
analysis  down  into  groups,  where  we  hold  some  parameters  constant,  we  begin  to  see  a 
pattern  emerge.  For  example,  from  comparison  of  heat-salt  (Figure  12)  with  cases  1  and 
2,  where  x  =  0.01  and  the  Prandtl  number  changes  from  7,  0.1,  to  100,  it  can  be  seen  that 
indeed  diameter  size  increases  with  a  decrease  in  Prandtl  number.  We  can  also  examine 
cases  where  the  Prandtl  number  remains  constant  and  the  diffusivity  ratio  varies,  by 
comparing  heat-salt  and  Case  3  (where  x  increases  from  0.01  to  0.1).  It  is  readily 
observed  that  increased  x  results  in  increased  finger  scales. 
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Temp  fiiihiif  (\\y) 


Pr  =  7.  x  =  0.1 


Case  4 

i 

■  “ 

'  ! 
ifl 

1  , 

I  V 

1  J 

1  i  < 

Pr=  1000.x  =  1/3 


Figure  18.  DNS  model  of  vertical  temperature  (fastest  rate  of  diffusion)  in  (x,  z) 

coordinates. 
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In  addition  to  what  has  been  observed  in  the  temperature  field  (horizontal)  section 
example  (Figure  16),  we  can  now  examine  finger  length  scales  and  their  dependencies 
from  the  vertical  sections  of  the  temperature  fields  (Figure  18).  It  is  clear  that  thinner, 
longer  fingers  occur  for  higher  Prandtl  numbers.  This  graphic  shows  that  a  decrease  in 
Prandtl  number  and  an  increase  in  the  diffusivity  ratio  both  cause  a  less  defined  finger 
formation. 
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Salinity  (x.y) 


. 


Pr=7.x  =  0.1 


Figure  19.  DNS  model  of  vertical  salinity  (slowest  rate  of  diffusion)  in  (x,  z)  coordinates. 
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Visual  inspection  of  the  vertical  salinity  sections  in  Figure  19  supports  the 
foregoing  inferences.  Comparing  heat-salt  with  Cases  2  and  4  implies  that  an  increase  in 
Prandtl  number  results  in  thinner,  longer  individual  fingers.  Magma  chamber  fingers 
(with  large  Prandtl  values)  have  been  observed  in  laboratory  experiments  as  being  slim  or 
slender  (Tait  1989)  and  is  verified  by  these  results.  A  more  objective  measure  of  the 
previously  observed  patterns  of  salt  fingers  is  given  by  the  aspect  ratio: 


Ar  = 


<T2  +Ty  +T72  > 
3<Tz2> 


(24) 


This  quantity  reduces  to  one  in  the  limit  of  isotropic  fingers  and  increases  as  fingers 
become  more  elongated.  Comparison  of  heat  salt  and  Case  3  (one  order  of  magnitude 
increase  in  diffusivity  ratio)  shows  only  a  slight  increase  in  the  aspect  ratio  (Ar).  A 
summary  of  the  dependencies  of  AR  is  shown  in  Figure  20  and  in  Table  4. 


Figure  20.  Aspect  Ratio,  Equation  (24),  for  Rp  =  1.9  for  Heat  Salt  and  Cases  1  to  4. 
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Aspect 

Ratio  at  R 

p=1.9 

Heat  Salt 

Case  1 

Case  2 

Case  3 

Case  4 

1.4928 

1.4349 

6.5011 

1.7506 

12.5670 

Table  4.  Aspect  Ratio  of  Heat  -  Salt  and  Cases  1  through  4. 


For  the  same  density  ratio,  Rp  =  1.9,  (Figure  20  and  Table  4)  the  aspect  ratio 
increases  with  an  increase  in  Prandtl  number.  As  x  increases,  the  aspect  ratio  increases 
slightly  as  reflected  by  the  differences  between  the  heat  salt  case  and  Case  3. 


Figure  21.  Case  1  Aspect  Ratio  as  a  function  of  Density  Ratio. 


Examination  of  Case  l(Figure  21  and  Table  5),  in  which  Pr  and  x  are  held 
constant  throughout  a  wide  range  of  Rp,  reveals  that  an  increase  in  density  ratio  is 
accompanied  by  a  mono  tonic  increase  in  aspect  ratio. 
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Case  1 

RP  Ar 

1.1 

1.2261 

1.3 

1.3271 

1.5 

1.3453 

1.7 

1.4101 

1.9 

1.4349 

2.1 

1.4796 

2.3 

1.4489 

2.5 

1.5791 

2.7 

1.6038 

2.9 

1.6564 

Table  5.  Aspect  Ratio  as  a  function  of  density  ratio  values  in  Case  1. 
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B. 


TRANSPORT  CHARACTERISTICS 


Equilibrium  Transport  -  Case  1  -  Rhol.9 


Dm* 


Equilibrium  Transport  -  Case  2  -  Rhol.9 


Dme 


Equilibrium  Transport  -  Case  3  -  Rhol  .9 
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Equilibrium  Transport  -  Case  4  -  Rhol.9 


Figure  22.  Temperature  and  salinity  fluxes  as  a  function  of  time  for  cases  1  through  4. 
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Case  #1 

Rp 

ft 

Fs 

Buoyancy 

Ratio 

1.1 

-6.4812 

-14.4047 

7.9235 

0.449936 

1.3 

-4.1979 

-9.8905 

5.6926 

0.424438 

1.5 

-3.1962 

-7.4937 

4.2975 

0.426518 

1.7 

-2.581 

-6.0696 

3.4886 

0.425234 

1.9 

-1.914 

-4.6501 

2.7361 

0.411604 

2.1 

-1.6215 

-3.8882 

2.2667 

0.417031 

2.3 

-1.4437 

-3.4342 

1.9905 

0.420389 

2.5 

-1.2994 

-3.0476 

1.7482 

0.426368 

2.7 

-1.0685 

-2.5368 

1.4683 

0.4212 

2.9 

-0.9298 

-2.1956 

1.2658 

0.423483 

Case  #2 

Rp 

ft 

Fs 

Buoyancy 

Ratio 

1.1 

-1105.08 

-1478.96 

373.8807 

0.7472 

1.3 

-468.524 

-705.2 

236.6764 

0.664384 

1.5 

-315.245 

-499.912 

184.6674 

0.6306 

1.7 

-242.235 

-398.488 

156.2523 

0.607887 

1.9 

-220.185 

-364.844 

144.6589 

0.603505 

2.1 

-196.214 

-328.582 

132.3677 

0.597154 

2.3 

-169.363 

-282.033 

112.6701 

0.600507 

2.5 

-168.507 

-286.134 

117.6272 

0.588908 

2.7 

-144.364 

-245.818 

101.4536 

0.587281 

2.9 

-149.324 

-250.356 

101.0315 

0.596448 

Case  #4 

Rp 

Ft 

Fs 

Buoyancy 

Ratio 

1.1 

-3009.76 

-3407.43 

397.6623 

0.883295 

1.3 

-1112.64 

-1321.63 

208.9909 

0.841869 

1.5 

-873.25 

-1047.29 

174.0426 

0.833817 

1.7 

-870.275 

-1032.55 

162.271 

0.842844 

1.9 

-769.48 

-885.398 

115.9173 

0.869079 

2.1 

-868.552 

-983.512 

114.9605 

0.883112 

2.3 

-650.989 

-720.552 

69.5631 

0.903459 

2.5 

-578.38 

-626.742 

48.3624 

0.922835 

2.7 

-328.609 

-344.985 

16.3752 

0.952534 

2.9 

-14.0786 

-14.4201 

0.3415 

0.976318 

Case  #3 

Rp 

Fj 

Fs 

Buoyancy 

Ratio 

1.1 

-234.914 

-324.005 

89.0909 

0.725032 

1.3 

-104.725 

-150.178 

45.453 

0.697339 

1.5 

-66.5649 

-98.8597 

32.2948 

0.673327 

1.7 

-49.8729 

-75.2535 

25.3806 

0.662732 

1.9 

-39.3424 

-60.2595 

20.9171 

0.652883 

2.1 

-33.3554 

-51.4704 

18.115 

0.64805 

2.3 

-26.356 

-41.2481 

14.8921 

0.638963 

2.5 

-12.2223 

-20.0053 

7.783 

0.610953 

2.7 

-15.3988 

-25.0506 

9.6518 

0.614708 

2.9 

-12.0376 

-19.733 

7.6954 

0.610024 

Table  6.  Tables  of  DNS  results  for  all  cases  including  average  temperatures  (FT) 
and  salinity  (Fs)  fluxes,  buoyancy  (Fs  -  FT)  and  flux  ratio  (FT/Fs). 


Typical  evolutionary  patterns  of  temperature  and  salinity  fluxes  (not  necessarily 
representing  temperature  and  salinity)  are  presented  in  Figure  22,  for  Rp  =  1.9.  The 
summary  of  all  calculations  in  Table  6  indicate  that  fluxes  are  most  strongly  affected  by 
variation  in  the  density  ratio  (Rp).  Our  limited  numerical  exploration  of  the  effects  of  the 
molecular  parameters  suggests  that,  for  any  given  Rp,  fluxes  increase  with  the  Prandtl 
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number  (Pr)  and  very  mildly  decrease  with  the  diffusivity  ratio  (t) — see  Tables  3  and  6. 
These  dependencies  are  consistent  with  the  prediction  based  on  the  growth  rate  balance 
model. 

The  purpose  of  the  following  calculation  is  to  detennine  whether  the  empirical 
parameterization  (Equation  (23))  developed  for  the  oceanographic  heat  salt  case  applies 
equally  well  to  other  configurations.  This  parameterization  was  applied  to  each  case, 
using  the  cftool  in  MATLAB.  Results  are  summarized  in  Table  7  (with  quality  of  the  fit 
indicated  in  the  last  column)  and  Figure  23. 


Empirical  Parameterization 

as 

bs 

aY 

bv 

cv 

Goodness 
Of  Fit 

Heat  Salt 

135.7 

-62.75 

2.708 

2.513 

0.5128 

99.7% 

Case  1 

5.113 

-0.7788 

427.7 

8.738 

0.4212 

95.9% 

Case  2 

503.9 

-160.6 

8.649 

3.663 

0.5928 

99.0% 

Case  3 

124.6 

-72.75 

0.359 

0.8303 

0.5772 

99.8% 

Case  4 

1189 

-482.2 

N/A 

N/A 

N/A 

93.5% 

Table  7.  Empirical  parameterization  of  the  temperature  and  salinity  fluxes. 
Parameterization  of  the  flux  ratio  for  Case  4  has  not  been  attempted  due  to  the 
nonmonotonic  pattern  of  y(Rp),  which  is  inconsistent  with  the  assumed  expression, 

Equation  (23). 


Results  in  all  cases  are  encouraging,  suggesting  the  universal  character  of  the 
assumed  parameterization  similar  to  Equation  (23),  as  quantified  by  the  goodness  of  fit 
for  all  simulations.  Examining  Table  7  shows  the  Prandtl  number  value  has  a  significant 
effect  on  the  magnitude  of  the  as  and  bs  coefficients:  as  Prandtl  number  increases  so  do 
the  magnitude  of  these  values. 
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Figure  23.  Best  fit  plots  for  Table  7,  Cases  1  through  4,  respectively.  Compare  to  Fleat 

Salt  Best  Fit  (Figure  15). 
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Buoyancy  Flux 


HeatSalt 
Casel 
Case  2 
Case  6 
Case  4 


R 


p 


Figure  24.  Buoyancy  flux  (Fs  -  Ft)  is  plotted  showing  that  as  Rp  increases  buoyancy 
decreases.  This  is  an  indirectly  proportional  relationship. 


One  of  the  key  characteristics  of  fully  developed  salt  fingering  is  related  to  its 
ability  to  transport  buoyancy  in  the  vertical  direction.  This  characteristic  is  thought  to  be 
involved  in  the  generation  and  maintenance  of  thermohaline  staircases  (Merryfield, 
2000).  The  buoyancy  flux  values  of  the  cases  presented  in  this  study  (Fs  -  FT)  are 
examined  in  Figure  24,  which  reveals  that  an  increase  in  Rp  results  in  a  decrease  in  the 
buoyancy. 
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Buoyancy  Flux 


Figure  25.  Buoyancy  flux  for  each  case  at  Rp  =  1.9. 


By  taking  a  closer  look  at  all  cases  at  a  fixed  density  ratio,  Rp  =  1.9,  (Figure  25 
and  Table  8)  we  can  see  the  two  largest  Prandtl  number  values  in  cases  2  and  4  (Pr  =100 
and  Pr  =  1000  respectively)  have  the  largest  buoyancy  values.  In  comparing  the  heat  salt 
and  Case  3  in  which  only  the  diffusivity  ratio  changes  (x  =  0.01  and  x  =  0.1,  respectively) 
there  is  a  slightly  smaller  buoyancy  value  for  Case  3.  Finally,  in  comparison  of  heat  salt 
and  Case  1,  with  a  decrease  in  Prandtl  number,  the  buoyancy  decreases. 


Buoyancy  flux  at  Rn=1.9 

Heat  Salt 

Case  1 

Case  2 

Case  3 

Case  4 

35.0034 

2.7361 

144.6589 

20.9171 

115.9173 

Table  8.  Buoyancy  flux  at  Rp  =  1.9  of  all  cases. 
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Figure  26.  Flux  ratio  in  all  case  studies. 


Another  critical  parameter,  also  implicated  in  the  formation  of  thermohaline 


staircases,  is  given  by  the  flux  ratio 


Y  =  - 


.  The  flux  ratio  generally  decreases  with  an 


Ls  y 


increase  in  Rp  as  seen  in  Figure  26.  However,  the  two  cases  (3  and  4)  with  larger 
diffusivity  ratio  (t)  show  an  increase  in  the  flux  ratio  for  higher  values  of  density  ratio 
(Rp).  Figure  27  and  Table  9  represent  the  flux  ratio  for  each  case  at  a  fixed  value  Rp  = 
1 .9,  which  also  indicate  that  the  largest  values  of  flux  ratios  occur  for  the  largest  value  of 
diffusivity  ratio. 


51 


Flux  Ratio 


Figure  27.  Flux  Ratio  for  all  cases  at  Rp  =  1.9. 


Flux  Ratio  at  Rn=1.9 

Heat  Salt 

Case  1 

Case  2 

Case  3 

Case  4 

0.538696 

0.4116 

0.603505 

0.65288 

0.8690788 

Table  9.  Flux  Ratio  for  all  cases  at  Rp  =  1.9. 


C.  COMPARISON  WITH  THE  THEORETICAL  MODEL 

The  relevant  values  of  the  adjustable  coefficient  (Q  vary  only  weakly  between 
various  cases  considered  in  this  study  (see  Figures  28  through  30).  The  lowest  value  of 
0=1.6  (obtained  for  Case  1)  is  not  too  different  from  that  largest  value  of  C=2.7  (heat-salt 
case)  considering  the  range  of  governing  parameters  and  fluxes.  Thus,  numerical 
evidence  in  this  study  supports  the  growth  rate  balance  theory  as  a  conceptual  model  of 
the  equilibrium  salt  fingering,  capable  of  predicting  dependencies  of  the  equilibrium 
temperature  /  salinity  transports  within  a  factor  of  two. 
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Figure  28.  Case  1 :  Temperature  (left  panel)  and  salinity  (right  panel)  fluxes  as  a  function 
of  density  ratio  (Pr=  0.1,  x  =  0.01).  Solid  curves  indicate  prediction  of  the 
theoretical  three-dimensional  model  and  DNS  values  are  represented  by  plus  signs. 


Figure  29.  Case  2:  Temperature  (left  panel)  and  salinity  (right  panel)  fluxes  as  a  function 
of  density  ratio  (Pr=  100,  t  =  0.01).  Solid  curves  indicate  prediction  of  the 
theoretical  three-dimensional  model  and  DNS  values  are  represented  by  plus  signs. 
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Figure  30.  Case  3:  Temperature  (left  panel)  and  salinity  (right  panel)  fluxes  as  a  function 
of  density  ratio  (Pr=  7,  x  =  0. 1).  Solid  curves  indicate  prediction  of  the  theoretical 
three-dimensional  model  and  DNS  values  are  represented  by  plus  signs. 
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VI.  DISCUSSION  AND  CONCLUSIONS 


Double-diffusion  is  an  essential  mixing  process  that  may  very  well  drive  the 
larger  scale  oceanographic  circulations  and  a  more  thorough  study  of  this  phenomenon 
may  lead  to  higher  accuracy  in  climate  and  ice  melt  models  and  global  and  regional 
oceanographic  forecast  models.  Understanding  this  process  requires  a  series  of  direct 
numerical  simulations  that  would  not  have  been  possible  without  parallel  processing. 
The  processing  time  utilized  in  our  study  totals  over  800,000  hours — 350,000  hrs  on 
Department  of  Defense  (DoD)  Systems  and  the  remaining  on  NPS  computers.  In  other 
words,  it  would  take  almost  a  full  century  to  complete  these  calculations  on  a  single 
processor.  The  parallelization  of  our  de-aliased  pseudo-spectral  model  allowed  us  to 
utilize  a  message  passing  interface  to  effectively  resolve  the  complex  calculations  of 
fluxes,  and  therefore  derive  an  understanding  of  how  these  micro-scale  mixing  processes 
occur  in  nature.  Such  an  understanding  will  facilitate  our  future  ability  to  model  these 
interactions  and  their  effects  on  large  scale  circulations.  This  too  requires  increased 
speeds  of  computer  processors  in  order  to  process  even  more  data  within  a  timely, 
forecastable  time  period. 

Our  model  predicts  the  following.  Double-diffusive  transport  of  heat  and  salt 
rapidly  intensifies  with  decreasing  density  ratio.  Fluxes  are  less  sensitive  to  molecular 
characteristics,  mildly  increasing  with  Prandtl  number  (Pr)  and  decreasing  with 
diffusivity  ratio  (x).  The  aspect  ratio  of  individual  fingers  monotonically  increases  with 
Prandtl  number,  with  tall  slender  fingers  realized  for  large  Pr  and  irregular  chaotic  fingers 
for  low  Pr.  Numerically  deduced  dependency  of  temperature  and  salinity  fluxes  on 
density  ratios  are  adequately  described  by  empirical  Equation  (23),  with  goodness  of  fit 
greater  than  93%  in  all  calculations.  Flux  ratio  shows  an  increase  at  larger  Rp  values  for  x 
>0.1  while  constantly  decreasing  for  x  =  0.01. 

In  addition  to  the  insights  based  on  the  numerical  modeling,  this  study  presents  a 
phenomenological  theory  for  the  vertical  transport  by  salt  fingers  in  unbounded 
temperature  and  salinity  gradients,  which  is  based  on  the  growth  rate  balance  Equation 
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(1).  By  utilizing  DNS  we  were  able  to  plot  equilibrium  fluxes  of  temperature  and  salinity 
as  a  function  of  time  and  explore  their  dependency  on  growth  rates.  This  dependency  led 
us  to  the  discovery  of  the  adjustable  coefficient  (C)  that  varies  slightly  for  a  wide  range  of 
parameters  and  fluxes.  Theory  assumes  that  the  statistical  equilibrium  is  reached  when 
the  growth  rate  of  salt  fingers  becomes  comparable  to  the  growth  rate  of  their  secondary 
instabilities.  This  assumption  can  be  rationalized  as  follows.  If  the  growth  rate  of 
primary  salt  fingers  exceeds  that  of  secondary  instabilities,  then  the  primary  instabilities 
are  unlikely  to  be  affected  by  perturbations  and  continue  to  grow.  If,  on  the  other  hand, 
the  secondary  instabilities  grow  much  faster  than  fingers,  they  would  rapidly  disrupt  the 
primary  modes,  dramatically  reducing  their  amplitude  and  transport  characteristics. 
Thus,  one  is  led  to  a  conclusion  that  the  approximate  balance  between  the  growth  rates  of 
primary  and  secondary  instabilities  is  inevitable  for  statistical  equilibrium  and  should  be 
satisfied  for  a  wide  range  of  governing  parameters. 

It  should  be  understood  that  the  growth  rates,  both  primary  (ki)  and  secondary 
(^2),  can  only  be  estimated  in  the  order  of  magnitude  sense.  Calculation  of  k|  is  made  for 
the  fastest  growing  salt  fingers  -  the  vertically  uniform  elevator  modes  -  whereas  fully 
developed  double-diffusive  convection  is  generally  dominated  by  irregular  structures  of 
finite  vertical  extent.  Similar  problems  arise  for  secondary  instabilities:  k2  is  estimated 
for  the  steady  and  vertical  background  field,  neither  of  which  is  an  exact  statement.  In 
view  of  these  difficulties,  we  only  claim  that  the  growth  rates  estimated  by  our  theory 
should  be  comparable,  but  not  necessarily  equal.  Thus,  we  consider  a  weak  form  of  the 
growth  rate  balance  in  Equation  (1).  The  nondimensional  order  one  coefficient  C  in 
Equation  (1)  cannot  be  determined  internally  on  the  basis  of  theory  and  therefore  it  was 
calibrated  using  direct  numerical  simulations  covering  a  wide  range  of  governing 
parameters.  For  each  (Pr,x),  theory  adequately  describes  the  (rapidly  decreasing) 
dependence  of  fluxes  of  diffusing  components  on  density  ratio.  The  calibrated 
coefficients  C  are  fairly  stable.  In  all  three-dimensional  simulations  -  simulations  that 
cover  three  orders  of  magnitude  in  Pr  and  one  order  in  x  -  the  values  of  C  are  limited  to  a 
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relatively  narrow  range  1.6  <  C  <  2.7.  The  theoretical  prediction  for  C  =  2.2  can  explain 
all  the  numerical  fluxes  in  this  paper  within  a  factor  of  two.  This  can  be  readily 
rationalized  using  the  growth  rate  balance  Equation  (1). 

Finally,  it  is  important  to  emphasize  the  theoretical  significance  of  our  model. 
We  believe  that  this  study  opens  opportunities  for  purely  analytical  explorations  that  are 
focused  on  the  growth  rate  balance  Equation  (1),  rather  than  on  the  more  complicated 
original  Navier-Stokes  system.  Our  optimism  is  partially  based  on  the  pronounced  low- 
order  dynamics  of  our  model.  Thus,  it  is  likely  that  the  general  character  of  our  solutions 
could  be  captured  by  an  analytical  low-dimensional  model. 
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VII.  RECOMMENDATIONS  FOR  FUTURE  WORK 


Imagine  our  children  and  grandchildren  planning  experiments  on  Jupiter  or  other 
celestial  bodies.  The  work  we  have  done  here  could  play  a  pivotal  role  in  how  they 
design  their  sensors  and  plan  their  missions.  As  we  have  no  easy  way  to  achieve  direct 
measurements  of  these  planetary  systems,  reliable  data  is  often  scarce  and  the  spatial 
location  of  the  elements  that  need  to  be  measured  must  be  pinpointed  to  the  most 
scientifically  relevant  areas.  Environmental  characteristics  of  any  future  planetary 
system,  need  to  be  pinpointed  to  allow  the  maximum  effectiveness  of  these  future 
missions. 

Even  closer  to  home  imagine  forecast  models,  atmospheric  and  oceanographic, 
with  successful  output  out  to  30  days  and  beyond.  Use  of  a  simple,  computationally 
compact,  algorithm  for  parameterization  of  micro-scale  mixing  in  general  and  salt  fingers 
in  particular  could  be  implemented  in  future  forecast  models  would  save  both  processing 
time  and  money.  Further  exploitation  of  this  microscale  feature  and  its  impact  on  large 
scale  circulations  can  greatly  improve  climate  and  operational  oceanographic  models.  Ice 
melt  and  climate  models  could  also  gain  accuracy  with  the  implementation  of  double 
diffusive  convection.  As  computer  technology  improves  and  processing  speeds  increase 
this  less  computationally  taxing  algorithm  can  easily  be  input  into  these  forecast  models. 

Successfully  modeling  double  diffusive  convection  could  impact  models  on 
magma  chamber  flows  and  may  even  prove  vital  in  volcanic  eruptions  and  possible 
earthquake  predictions.  With  successful  earthquake  and  volcanic  eruption  predictions  we 
can  also  predict  tsunamis  and  have  more  warning  time  in  order  to  prepare  for  these  often 
catastrophic  events.  The  events  in  the  past  few  years  show  that  these  advanced  warning 
systems  are  needed  to  prevent  the  mass  loss  of  life  with  the  “Ring  of  Fire”  surrounding 
the  Pacific  Ocean. 
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